Time-dependent, four-point density correlation function description of dynamical 
heterogeneity and decoupling in supercooled liquids 
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Dynamical heterogeneity and the decoupling of diffusion and relaxation in a supercooled liquid 
is investigated via a time-dependent, four-point density correlation function. We show that the 
main contribution to the corresponding generalized susceptibility X4(*) in molecular dynamics 
simulation of a Lennard- Jones liquid arises from spatial correlations between temporarily local- 
ized ("caged") particles. By comparing X4(^) with a generalized susceptibility XM{i) related to a 
correlation function for the squared particle displacements, we demonstrate a connection between 
dynamical heterogeneity and the decoupling of relaxation and diffusion. 
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Spatially heterogeneous dynamics ( "dynamical hetero- 
geneity") in otherwise homogeneous supercooled, gl^^ 
forming liquids is now well established in experimentsETEl 
close to the glass transition temperature Tg, and this het- 
ergeneity is even apparent at higher temperatures a^*?*^ 
the mode couplingl^ temposatuee Tc in simulationsLrtEl. 
For example, recent studiesETtjO of the dynamics of su- 
percooled, glass-forming polymeric and binary simple liq- 
uids in terms of the correlations of monomer or particle 
displacements revealed the dynamical heterogeneity of 
these liquids and a rapidly growing range of correlated 
motion on cooling towards Tc- At the same time, the de- 
coupling of translational diffusion and relaxation as well 
as translational and rotational diffusion in [these fluids 
is also well known, and simulations showEJllj that this 
decoupling begins well above Tc where dynamical hetero- 
geneity first appears. Several authors have argued that 
the decoupling of diffusion and relaxation is a direct re- 
sult of dynamical heterogeneity, with the slowest particles 
dominating structural -cela^alinu and the fastest particles 
dominating diffusionldfflElmO. 

In this Letter we use a four-point time correlation 
function of the density to probe dynamical heterogene- 
ity in a glass-forming liquid, and elucidate the connec- 
tion between this heterogeneity and the decoupling of 
bulk transport processes. This four-point function was 
first iavestigated in a supercooled-Jjquid by Dasgupta, 
et al.Ea, and recently Donati, et al.Ej have demonstrated 
analytically and computationally the interesting behavior 
of the related generalized four-point-|Susceptibility Xiit) 
(defined below). As shown in Ref.ll3, Xii^) can be rep- 
resented in terms of the fluctuations of an "order pa- 
rameter" that is a bilinear, time-dependent product of 
densities. Here we show that the self-part of X4(t) is di- 
rectly related to spatial correlations between temporar- 
ily localized particles, while the distinct part is related 
to the correlated motion of particles into positions pre- 
viously occupied by neighboring particles. We evaluate 
these quantities for a cold Lennard- Jones (LJ) liquid, and 



show that in this system Xii^) is dominated by growing 
spatial correlations between localized particles. We then 
compare the behavior of Xii^) with a generalized time- 
dependent susceptibility related to a correlation function 
of squared particle displacements. From these two quan- 
tities we find two different characteristic time scales: the 
time scale on which temporarily localized particles are 
most spatially correlated scales with temperature like the 
structural relaxation time, while the time scale on which 
the correlation between squared particle displacements is 
strongest scales like the inverse diffusion coefficient. In 
this way, we demonstrate that the decoupling of diffusion 
and relaxation in this model liquid arises from dynamical 
heterogeneity. 

Consider a liquid of N particles in a volume V , with 
density p{r,t) = J2iLi ~ ^ii^))- ^he simplest density 
correlation function that contains information on corre- 
lated particle motion is fourth order. We write this func- 
tion in terms of the deviations of p(r, t) from its average 
value, Ap{r,t) = p{r,t) - po, where po = (p) = N/V, 
and (...) denotes an ensemble average: 

.F4(ri,r2,t) = (Ap(ri,0)Ap(ri,i)Ap(r2,0)Ap(r2,t)) 

-(Ap(ri, 0)Ap(ri, t))(Ap(r2, 0)Ap(r2, i)). (1) 

Terms involving one position only are subtracted in 
Eq. (0) since they contain no information on spatial cor- 
relations of particle motions. ^4(ri,r2,t) can be written, 

:?^4(ri,r2,<) = tJ4(ri,r2,<) + AJ^4(ri, r2, i), 

where the two-point, two-time, fouijth.-nrder correlation 
function of densities Qi is defined asEflEfl 

^4(ri,r2,t) = (p(ri,0)p(ri,<)p(r2,0)p(r2,i)) 

-(p(ri,0)p(ri,i))(p(r2,0)p(r2,0). 

AjF4(ri, r2, t) consists of second- and third-order correla- 
tion functions of density. A straightforward calculation 
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shows that J J dridr2AT4{ri,r2,t) vanishes by symme- 
try, and as a result, the volume integrals of ^4(ri,r2,t) 
and S4(ri,r2,t) are equal to each other and correspond 
to the same generalized susceptibility X4(^)i 

It is straightforward to show that X4 {t) can be written as 



xl{t)^^[{Ql{t))-{Qo{t)?], 



(2) 



where [3 — l/fc^T, and the time-dependent "order pa- 
rameter" Qo(i) equals 



N N 



Qo{t)= /'rfrp(r,0)p(r,t)-^^<5(r,(0)-r,(0). (3) 

i=i 3=1 

In a simulation, QQ{t) is numerically ill-defined (for a 
finite system) since the probability that particle j exactly 
replaces particle i is infinitely small. Following ParisiE2l, 
we therefore modify Qo(t) by an "overlap" function w(r) 
that is unity inside a region of size a and zero otherjSSiise, 
where a is taken on the order of a particle diameterEBEil. 
This leads to an a-dependent counterpart to Qo(t), 

Q{t) = j rfridr2p(ri, 0)p(r2, i)w(|ri - r2|). 

N N 

= EE / dvw{\v\)5{v + vm-^At)) 

A — 1 „■ — 1 ^ 



i=l 3 = 1 • 
N N 

i=i 3=1 



(4) 



where = (0) — rj (0) and fli = (t) — (0) is the 
displacement of particle i during the time [interval from 
zero to t. We choose a — O.SaAA as in Ref.E3. 
Replacing Qo{t) in Eq. (p|) by Q{t) yields 



(5) 



which gives the following expressionlla for Xii^) terms 
of the four-point correlation function G4(ri, r2, ra, r4, i): 



I3V f 

^4(0 = ^ dridr2dr3dr4w{\ri - rsD^drg - r4| 



X G4(ri,r2,r3,r4,t) 



(6) 



where 



G'4(ri,r2,r3,r4,i) = (p(ri, 0)p(r2, i)p(r3, 0)p(r4, t)) 

-(p(ri,0)p(r2,t))(p(r3,0)p(r4,t)). (7) 

We can write Q in terms of its self and distinct parts, 
Q — Qs + Qd- The self part Qs corresponds to terms 
with i = j in Eq. (Q): 
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(8) 



where fXj is the magnitude of jlj. The distinct part Qd 
is equal to 



N N 



(9) 
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Then X4(i) can be decomposed into self (xss), dis- 
tinct ixDo), and interference (xsd) parts: x = XSS + 
Xdd + XSD- From Eq. (5), xss and xdd describe the 
fluctuations of Qs and Qd, respectively, and xsd de- 
scribes the cross fluctuations: xss {Q%) ~ {QsYi 
Xdd oc - (<3^^' and xsd oc (QsQd) - {Qs){Qd)- 
According to Eq. (||), Qsit) contains only contributions 
from small displacements, fii < a, since w{^i) = for 
fi'i > a, and thus xss{t) is the susceptibility of localized 
particles, those which during a time interval [0,t] move 
less than a distance a. In contrast, Qoit) contains con- 
tributions from particles for which \rij — flj\ < a; that 
is, particles that replace neighboring particles in a time 
interval [Qjt]. 

In Ref.B, a different generalized susceptibility Xu{t) 
was defined in terms of the fluctuations in an "order pa- 
rameter" given by the total particle displacement in a 
time interval t, U{t) = J2iLi — J dru{r,t), where 
the displacement density field u{r,t) — J2i=i ~ 
rj(0)). In a similar fashion, Eq. (||) can be rewritten as 
Qs{t) — J drqs{r,t) where the "localization density" 
field qs{r, t) is qsir, t) = J2iLi w{fii{t))6{r - ri(0)). Here 
we compare Xi{t) with XM{t), defined as 



XM{t) = 



{M{t)y 



where M{t) = J2tLi t^Ut) (i-e- M{t) is the sum of the 
squared displacements for one system in a time interval 
[0,t]). Like Xu{t)j XM{t) is proportional to the volume 
integral of a correlation [fanction of (in this case squared) 
particle displacementsQa. Note that both the displace- 
ment density field u{r, t) and squared-displacement den- 
sity field m(r, t) are dominated by delocalized, or mobile 
particles, while the localization density field qs{j^,t) is 
dominated by localized, or immobile particles. 

To evaluate these quantities we use data obtained from 
a molecular dynamics simulation of a model LJ glass- 
former. The system is a three-dimensional binary mix- 
ture (50:50) of |5Q0 particles interacting via LJ interac- 
tion parametersES. We analyze data from state points at 
seven diffdjent temperatures T approaching Tc ~ 0.592 
from aboveE^ at a constant density p w 1.3. (In the re- 
mairuder of this paper, all values are quoted in reduced 
unitsc3.) All quantities presented here are evaluated in 
the NVE ensemble following equilibration of the system 
at each gtate point. Further simulation details may be 
found inN-El. 



■[{M^{t))-{M{t)Y^ 



(10) 
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In Fig. |l|a the susceptibility X4(i) calculated via Eq. ^ 
is shown as a function of time for differentr-walues of T. 
As found for a different LJ mixture in Ref.ll3, for all T, 
Xiit) is zero at short time and attains a small constant 
value at large time, and has a maximum at some inter- 
mediate time t|. Both t| and the amplitude of the peak, 
XiitDi increase strongly with decreasing T. At the low- 
est value of T, the amplitude of Xii^) decreases, possitfe 
due to finite size effects or to the change in dynamicsEj 
near T^. The inset shows the self, distinct and cross terms 
of Xiit) for one value of T, and we see that xss is in- 
deed the dominant term. Thus, x^it) is dominated by 
the growing range of spatial correlations between local- 
ized particles in this fluid, and t| is the time when this 
correlation is strongest. In fact, several authors have re- 
ported evidence of a growing length scale associated with 
solid-like behavior in dense fluida^J. 




FIG. 1. (a) Time dependence of the susceptibility Xii't) a-t 
various temperatures, as indicated in (b). Inset: Self, distinct 
and cross terms of X4(^) at T = 0.62. (b) Time dependence of 
the "squared-displacement" susceptibility XM{t) at the same 
values of T as in (a) . 

Fig. lb shows XM{t) calculated from Eq. |l^, as a func- 
tion of time for different values of T. We find that 
XM{t) becomes negligable at small and large times and 
has a maximum at some intermediate time where the 
spatial correlation of squared particle displacements is 
strongest. This behavioc is similar to that exhibited by 
Xu(t) calculated in Ref.S. 

As shown in Fig. 2, both Xii^l) and XmI^m) increase 
strongly with decreasing T (with the exception of XiitX) 
at the lowest temperature). Over the limited temper- 
ature range of our simulations, both functions may be 



reasonably fitted by power-law functions {T — Tc)~^ with 
Tc = 0.592, with the apparent exponents 74 = 0.80±0.07 
and 7A/ = 0.87 ± 0.05, as shown in the figure. (In fitting 



to the value Tc = 0.592 
) Of course, precise de- 



the power law, Tc is held fixi 
determined in previous world 
termination of the functional form requires simulations 
at lower temperatures, and larger simulations to reduce 
any possible finite size effects expected due to the grow- 
ing range of correlated particle motion and localization 
driving the growth of XM{t) and X4(t), respectively. 
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FIG. 2. Temperature dependence of x^i^-t) ^nd XM{t*M). 
The solid and dashed lines are power law fits to the data 
as indicated (excluding the lowest temperature). The error 
bars are estimated from deviations between three indepen- 
dent samples, where for each sample, X4(^) a-^d XM{t) are 
calculated by averaging over 128 independent time origins. 

Figs. 3(a) and (b) show the T-dependence of tl and 
t*j^j, respectively, and compare them with both the in- 
verse self-diffusion coefficient D^^ and the structural re- 
laxation time Tq. (Here D is calculated from the mean 
square displacement for the B (small) particles, and Ta 
is calculated by fitting the a-relaxation part of the self 
intermediate scattering function for the A (large) parti- 
cles (not shown, see Ref.t^l) by a stretched exponential 
function.) Also shown are power law fits to^D^^ and 
Ta excluding the lowest temperature (see Ref£j for de- 
tails). Diffusion and relaxation are found to be "decou- 
pled" in this cold liquid, aflr^|)tfP5¥|?f^ i^ many other real 
and simulated cold liquidsBcHlj^tj'EZl. In the present sys- 
tem, we find that with = 0.592, -fo = J-Jl ± 0.03 
and jr = 1-41 ± 0.07. Remarkably, we findEa that the 
T-dependence of i| coincides within our numerical er- 
ror with that of Ta, while t^^^ behaves like D^^. That 
is, the time scale on which the localized particles are 
most spatially correlated scales with temperature like the 
structural relaxation time, and the time scale on which 
the correlation between squared particle displacements 
is strongest scales like the inverse diffusion coefficient. 
Thus, our data demonstrates that the "decoupling" of 
diffusion and relaxation (or viscosity) may be directly 
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attributed to the emergence of dynamical heterogeneity, 
as arguedpby, e.g., Siilescu and cowMkeraJ, Ediger and 
coworkersQ, StillingeiM, and Douglasli3. 
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FIG. 3. (a) Temperature dependence of Tq and the time tl 
at which X4(i) exhibits a maximum, (b) Temperature depen- 
dence of the inverse self-diffusion coefBcient D and the time 
tlf at which XA/(i) exhibits a maximum. The solid lines are 
power law fits to Tq, and D^^ respectively (excluding the low- 
est temperature), with Tc fixed. Insets: Comparison of t1 and 
f^i with both and r^. As plotted, a line of zero slope 

(dashed line) indicates proportionality. 

Our results demonstrate the importance of time- 
dependent higher order density correlation functions in 
the characterization of dynamical heterogeneity in super- 
cooled liquids, and the ramifications of this heterogeneity 
for the bulk dynamics. In particular, the increasing am- 
plitude of the generalized time-dependent susceptibility 
Xiit) with decreasing T, as shown also in Ref.E£l, demon- 
strates an essential difference between two- and four- 
point density correlation functions in these fluids. For a 
glass-forming LJ liquid, we have shown that Xii^) is dom- 
inated by growing spatial correlations between temporar- 
ily localized particles. Finally, we have demonstrated 
that the decoupling of diffusion and structural relaxation 
observed in supercooled liquids follows naturally from dy- 
namical heterogeneity, as discussed by Siilescu, Ediger, 
Stillinger and Douglas: the time scale for spatial cor- 
relations of localized particles to develop governs struc- 
tural relaxation, while the (different) time scale for the 
development of spatial correlations of squared particle 
displacements, which is dominated by mobile particles, 
governs diffusion. We note that it should be possible to 



determine the 4-point functions studied here in colloidal 
suspensions using particle tracking methods. 

We thank C. Donati, J.F. Douglas, S. Franz, R.D. 
Mountain, G. Parisi, P.H. Poole, and F. Starr for valu- 
able comments on the manuscript. 



* On leave from Institute of Automation & Electrometry, 
Russian Academy of Sciences, Novosibirsk, 630090, Russia. 
^ H. Siilescu, J. Non-Cryst. Solids 243, 81 (1999), and refer- 
ences therein. 

^ R. Bohmer, Curr. Opin. Sol. State Mater. Sci. 3, 378 (1998) 
and references therein. 

^ See, e.g., K. Schmidt-Rohr and H.W. Spiess, Phys. Rev. 
Lett. 66, 3020 (1991); A. Heuer, et al., Phys. Rev. Lett. 
95, 2851 (1995); R. Richert, J. Non-Cryst. Solids, 172- 
174, 209 (1994); R. Bohmer, et al., J. Non-Cryst. Solids 
235-237, 1 (1998); F. R. Blackburn, et al., J. Non-Cryst. 
Solids 172-174, 256 (1994). 

F. Fujara, B. Geil, H. Siilescu, and G. Fleischer, Z. Phys. 
B 88, 195 (1992); I. Chang, et al., J. Non-Cryst. Solids 
172-174, 248 (1994); H. Siilescu, Phys. Rev. E 53, 2992 
(1996); I. Chang and H. Siilescu, J. Phys. Chem. B 101, 
8794 (1997). 

^ M.T. Cicerone, F.R. Blackburn and M.D. Ediger, Macro- 
molecules 28, 8224 (1995); M.T. Cicerone and M.D. Ediger, 
J. Chem. Phys. 104, 7210 (1996). 

^ L. Andreozzi, A. Di Schino, M. Giordano, and D. Leporini, 
Europhys. Lett. 38, 669 (1997). 

See, e.g., T. Muranaka and Y. Hiwatari, Phys. Rev. E 51, 

R2735 (1995); R.D. Mountain, J.Chem.Phys. 102, 5408, 

(1995); M. Hurley and P. HarroweU, Phys. Rev. E 52, 

1694 (1995); D. Perera and P. HarroweU. J. Non-Cryst. 

Sol. 235-237, 314 (1998); B. Doliwa and A. Heuer, Phys. 

Rev. Lett. 80, 4915 (1998); A. Onuki and Y. Yamamoto, 

J. Non-Cryst. Sol. 235-237, 34 (1998). 
* C. Donati, P.H. Poole, and S.C. Glotzer, Phys. Rev. Lett., 

82, 6064 (1999); C. Bennemann, C. Donati, J. Baschnagel 

and S. C. Glotzer, Nature 399, 246 (1999). 
^ S. C. Glotzer and C. Donati, J. Phys.: Cond. Matt. 11, 

A285 (1999); W. Kob, C. Donati, P.H. Poole, S.J. Plimpton 

and S.C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997). 
^° C. Donati, S. Glotzer, P. Poole, W. Kob, and S. Pfimpton, 

Phys. Rev. E 60, 3107 (1999). 
" W. Gotze, L. Sjogren, Rep. Progr. Phys. 55, 241 (1992). 

R. Yamamoto and A. Onuki, Phys. Rev. Lett. 81, 4915 

(1998); Phys. Rev. E 58, 3515 (1998). 

W. Kob and H.C. Andersen, Phys. Rev. E 51, 4626 (1995). 
C. Bennemann, W. Paul, K. Binder, and B. Diinweg, Phys. 
Rev. E 57, 843 (1998). 

S. Kammerer, W. Kob and R. Schilling, Phys. Rev. E 56, 
5450 (1997). 

F. H. Stillinger and J. A. Hodgdon, Phys. Rev. E 50, 2064 
(1994); Phys. Rev. E 53, 2995 (1996). 

J.F. Douglas and D. Leporini, J. Non-Crystalline Sol. 235- 
237, 137 (1998). 



4 



' C. Dasgupta, A.V. Indrani, S. Ramaswamy and M.K 
Phani, Europhys. Lett. 15 307 (1991). 



C. Donati, S. Franz, G. Parisi, and S.C. Glotzer, cond 



mat/9905433 , 

G. Parisi, J. Phys. A: Math. Gen. 30, L765-L770 (1997). 
M. Cardenas, S. Franz, and G. Parisi, J. Phys. A: Math. 
Gen. 31, L163 (1998); J. Chem. Phys. 110, 1726 (1999). 
D. Lancaster and G. Parisi, J. Phys. A: Math. Gen. 30, 
5911 (1997). 

The LJ interaction parameters are asB joAA = 5/6, oab ~ 
(pAA + (Jbb)/2, and eAA = <^ab = ess. The masses are 
given by niB/mA — 1/2. The length of the sample is L = 
7.28aAA and the potential was cut and shifted at 2.5cra,/3. 
All quantities are reported in reduced units: T in units 
of tAA, lengths in units of ctaa and time in units of r = 
(mscr^A/48e)^'^^. Adopting Argon parameters for the A 
particles gives r = 3 x lO"'^^ s. 
24 T.B. Schr0der and J.C. Dyre, J. Non-Cryst. Sohds 235- 
237, 331 (1998). 



T.B. Schr0der , S. Sastry, J.C. Dyre and S.C. Glotzer, cond- 



mat/9901271 



R.D. Mountain, ACS Symposium Series 676, 122 (1997); 
R. Ahluwalia and S.P. Das, Phys. Rev. E 57, 5771 (1998); 
G. Johnson, et al., Phys. Rev. E 57, 5707 (1998). 
P. AUegr ini, J.F. Douglas an d S.C. Glotzer, Phys. Rev. E, 
in press (cond-mat/990731J ). i-i 
We also calculated xu{t) following Refs.El; we find that al- 
though t* scales more like than like Tq, t^j scales more 
like than does tZ. 



5 



